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ABSTRACT 


The Naval Postgraduate School mixed layer model is augmented to include an 
entrainment zone with finite thickness. The role of entrainment shear production of 
turbulent kinetic energy is investigated by comparing model results that include the new 
entrainment zone with observations at Ocean Weather Station Papa in the North Pacific 
and with model predictions that do net include the new entrainment zone feature. Al¬ 
though it is not yet clear that annual-period model forecasts are improved significantly, 
it is shown that the entrainment zone processes play a significant role in vertical fluxes 
and in the turbulent kinetic energx budget of the upper ocean under warming conditions. 
Furthermore, it is found that the improved entrainment zone more accurately reproduces 
the temperature gradients of transient thermoclines observed at OWS Papa. 



Accession For 


NTIS GRA&I 
OTIC TAB 
Unannounced 
Justification. 




□ 

□ 


H?trlbutlon/_ 

| Availabili ty Codes 
(Avail and/or" 
Special 





TABLE OF CONTENTS 


I. INTRODUCTION.1 

A. PURPOSE .1 

B. BACKGROUND .1 

C. APPROACH .4 

II. THEORY . 5 

A. GOVERNING EQUATIONS .5 

B. MODEL EQUATIONS . 6 


III. METHOD .. II 

A. RESPONSE TIME ANALYSIS .11 

B. PRELIMINARY ADJUSTMENTS .14 

C. TUNING THE MODEL .17 

IV. RESULTS .21 

A. OVERVIEW ..21 

B. SST SIMULATIONS .21 

C. UPPER OCEAN THERMAL STRUCTURE SIMULATION .24 

D. EXAMINING ENTRAINMENT SHEAR PRODUCTION .30 

V. SUMMARY AND REMARKS .43 

REFERENCES .46 


INITIAL DISTRIBUTION LIST 


48 























LIST OF TABLES 


Table 1. ANNUAL SST ERRORS. 

Table 2. EFFECT OF ENTRAINMENT SHEAR PRODUCTION 


v 






LIST OF FIGURES 


Figure 1. The NPS Model .2 

Figure 2. Numerical solution with N'iiler model .12 

Figure 3. Error arrays . 20 

Figure 4. Model Estimated SST .. 25 

Figure 5. Temperature profiles .28 

Figure 6. Temperature profiles ..29 

Figure 7. Temperature profiles using identical model constants .32 

Figure 8. Wind Stress . 39 

Figure 9. Heat Flux .. 40 

Figure 10. Difference in SST . 41 

Figure 11. Thermoclinc representation .42 


VI 













ACKNOWLEDGEMENTS 


The author wishes to express sincere appreciation to all who have made possible 
the completion of this research. 

Special thanks are extended to Dr. Bill Garwood for the time and effort he has 
spent assisting in the direction and analysis of this investigation, in providing the NPS 
model, in touble-shooting problems encountered during the research, and in preparing 
this manuscript. 

Gratitude is also extended to the following individuals for their generous assistance: 
Dr. Pecheng Chu, who provided encouragement, offered suggestions, and read the 
manuscript; Dr. Robert Haney, who assisted in interpreting some of the theory; Arlene 
Bird, who provided expertise in writing computer programs for plotting the data; and 
Robert Rcehm, who provided encouragement and typing services. 

Finally, the author would like to thank his beloved wife, Carla, whose continuous 
inspirational support and selfless sacrifices made the timely completion of this endeavor 
possible. 



I. INTRODUCTION 


A. PURPOSE 

The purpose of this study is to investigate possible improvements in the Naval 
Postgraduate School oceanic mixed layer model, which is used in prediction of the upper 
ocean thermal structure, by including a finite-thickness entrainment shear production 
zone at the base of the mixed layer and allowing entrainment mixing to occur below the 
well-mixed surface layer. 

B. BACKGROUND 

In this study the oceanic mixed layer is considered to oe a fully turbulent region 
bounded above by the air-sea interface and below by a dynamically stable water-mass. 
The surface layer is assumed to be homogeneous in temperature and salinity and to have 
nearly uniform horizontal currents, except for domains of large shear near the surface 
and in the entrainment zone. The upper shear zone is forced by surface buoyancy flux, 
wind-driven shear production of turbulence, and wind-driven wave action, while the 
lower is characterized by entrainment of non-turbulcnt dense water from below, This 
lower zone is the region of interest for our study. The general features of this three-layer 
model are presented in Figure 1 on page 2. 

As pointed out by Gaspar (1988), the vertical forcing of surface heat and momen¬ 
tum fluxes dominates the physics of the upper ocean. Consequently, the mixed layer can 
be reasonably treated as a bulk layer (often referred to as a slab) where all horizontal 
gradients are regarded as negligible. Kraus and Turner (1967) intioduced the first 
oceanic one-dimensional bulk model by vertically integrating the total turbulent kinetic 
energy equation over the depth of the layer. Other bulk models have since been intro¬ 
duced which are based on the same principal assumptions as the Kraus and Turner 
model, but they differ in their methods of parameterizing the physical processes 
associated with the mixed layer. The Naval Postgraduate School (XPS) model is a 
numerical solution to the Garwood (1977) model which was the first bulk model to 
formulate the turbulent kinetic energy budget into three separate components, using 
bulk second-order closure to solve for the vertical averages of these components. 
Reviews of various bulk models and their parameterizations are provided by 
Zilitinkcvich et al. (1979) and Garwood (1979). 
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Figure 1. The NFS Model: Treatment of physical processes as shown for a typ¬ 
ical temperature profile of the upper ocean. 


The capability of bulk models to simulate evolutions of the mixed layer has been 
demonstrated by several authors, including Martin (1985) and Caspar (1988). Although 
second order turbulence closure models that compute the vertical profiles of the turbu¬ 
lent kinetic energy have also been introduced (eg., Mcilor and Durbin 1975), the bulk 
models yield comparable results and are computationally more efficient. 

Like several other mixed layer models, earlier versions of the NTS (Garwood 1977) 
model have been evaluated against data collected at ocean weather stations. A review 
of several of these investigations is provided by Martin (1985). This study will use data 
gathered at Ocean Station Papa which was in operation in the eastern North Pacific 
between 1949 and 1981. 







Station Papa was occupied alternately by two vessels operated by the Marine 
Services Branch of the Canadian Ministry of Transport. Observations at the station 
consist of frequent subsurface bathythermograph (BT) casts and three-hourly surface 
meteorological measurements. The meteorological data include dry and wet-bulb tem¬ 
perature, sea surface temperature (SST), wind speed and direction, and fractional cloud 
cover in octals. The BT casts were conducted at intervals varying between twenty min¬ 
utes and several days and yielded vertical temperature profiles that were digitized at five 
meter intervals. All measurements were confined to the region between latitudes 49° N 
and 51° X and longitudes 144° W and 146° W. 

Detailed descriptions of oceanic and atmospheric conditions at Station Papa arc 
provided in Tabata (1961 and 1965). The horizontal currents are weak in this region. 
Additionally, ad\cction by upwelling is insignificant in comparison to horizontal 
advcction. These characteristics make this location ideally suited for testing one- 
dimensional mixed layer models. 

All previous tests of the NPS model at Papa have used numerical versions that 
neglected entrainment shear production and assumed the entrainment zone to be 
infinitesimally thin. These models did not allow the entrainment zone to have a thick¬ 
ness, <5, greater than the model grid size (Ar = 1 m). One of these versions of the model 
was evaluated by Martin (19S5) for annual simulations of mixed layer depth evolution 
and SST changes. Martin suggested that the model might be further improved for op¬ 
erational use if it included entrainment shear production along with the momentum 
budget ncssecary to include this process. 

De Szoekc and Rhines (1976) and Garwood and Yun (1979) suggested that in most 
circumstances the production of turbulent kinetic energy by entrainment shear pio- 
duction is negligible. Recent research, however, indicates that the process may actually 
be significant. Garwood (1987) has shown that to the first order, the entrainment zone 
thickness, <5, is independent of the surface buoyancy flux and is dependent largely on the 
Coriolis parameter and the surface wind stress. Furthermore, the intensity of turbulence 
and mixing in the entrainment zone is not dependent upon the Richardson number when 
the layer is at the critical value (7?/ f ,s:-j-). Rather, Garwood et al. (1989) have shown 

H 

that the intensity of mixing in the entrainment zone can be mostly controlled by the 
fluxes of energy from above. 
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C. APPROACH 

This study will investigate the significance of the entrainment shear production 
process in predictions of upper oceanic features, particularly sea surface temperature. 
Conclusions will primarily be sought by evaluating two versions of the N’PS model 
against Station Papa data. One version will include a formulation for the process, and 
the other version will neglect it. A comparison of the model results will then be pre¬ 
sented. Additionally, the time and depth dependence of the entrainment shear pro¬ 
duction process will be demonstrated theoretically by following the method of de Szocke 
and Rhines (1976). The evaluation procedure will be preceded by a review of the gov¬ 
erning equations and a presentation of the NPS model, highlighting features neglected 
in other studies. 
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II. THEORY 


A. GOVERNING EQUATIONS 

As presented in Garwood (1977), the NPS model is based on the Navier-Stokes 
equations of motion without geostrophic components, the first law of thermodynamics, 
the conservation of salt mass, the continuity equation for incompressible water, and an 
equation of state based on temperature and salinity. It assumes that the turbulent kinetic 
energy of the overlying mixed layer is transformed into potential energy as the lower 
water-mass is destabilized and entrained into the upper layer. Applying these assump¬ 
tions, the model uses exchanges of heat, salt, and momentum to estimate changes in the 
heat content and turbulence of the mixed layer. This estimate is then used to compute 
a change in temperature and a corresponding deepening or shallowing rate of the mixed 
layer. 

This presentation of the model will employ rectangular coordinate axes with .v 
positive to the east, y positive to the north, and z positive upward originating at the 
ocean surface. Likewise, the eastward, northward, and upward components of velocity 
will be represented as u , r, and w, respectively. 

Momentum is assumed to be imparted only by wind forcing. The existence of a 
horizontal pressure gradient is assumed to have no effect on vertical shear for creation 
of turbulent kinetic energy and is neglected. Consequently, wind stress components, 
r, and x y , determine the values of the horizontal friction velocity: 

«* = [( t » 2 + ( Vp ) 2 ]’ /4 ( 1 ) 

where p denotes the density of sea water. 

Surface salinity flux is due to net evaporation minus precipitation (e-p): 

IV(0 ) = -S(e-p) (2) 

where an overbar or an upper case letter denotes a vertical mean quantity and a prime 
represents a turbulent or fluctuating quantity. This notation follows the convention of 
representing a given variable * as: x = X+x' , where the mean across the mixed layer h 
is defined: 
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( 3 ) 


X = 



Heat fluxes are defined by the following equation: 

Qo~Qs~Qb~Qh~ Qe 

where Q e is net downward heat flux, Q s downward solar radiation, Q b upward back ra¬ 
diation, Q„ upward sensible heat flux, and Q E upward evaporative heat flux. A complete 
description of how these fluxes are determined for the model is provided by Garwood 
(1976), and an analysis of the parameterization of the absorption of solar radiation in 
the upper ocean is provided by Gallacher et al. (1983). 

The unique aspects of the N'PS model are manifested in its method of computing 
the rate of deepening or shallowing of the mixed layer. Of particular note, is the treat¬ 
ment of dissipation in a three-component sense. To introduce this feature of the model, 
it is best to begin with a brief disscussion of the model equations. A detailed description 
of the model is provided by Garwood (1976 and 1977). 

B. MODEL EQUATIONS 

Essential to the solution of this model is an accurate assessment of the entrainment 
velocity, w t . Letting: 


Hi 

ci 


= H' e - IV(-h). 


( 4 ) 


For Ocean Station Papa it is assumed that the mean vertical motion is negligible or 
1F( -/:) = 0. Thus the equation becomes: 


I!l = w 

£t e ' 

The entrainment velocity is found by solving 

rv T2\l 12 

E(\V ) ' 


\gh(ct&T - /IAS) + E\ 


( 5 ) 


( 6 ) 


where g is the acceleration of gravity, a and /? the thermal and saline expansion coefli- 
cients, and E= u' 2 + v' 2 + w' 2 . The variables AS and A7' denote the change in mean 
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salinity and temperature across the entrainment zone and are obtained using the jump 
condition. The principal assumption being applied is that the mixed layer is approxi¬ 
mately homogeneous in the mean properties of T and S. Often the effects of variations 
in heat and salinity are combined into a single buoyancy term, b. Thus the mean 
buoyancy jump across the entrainment zone 5 can be expressed as: 

A B = B{ - h ) - B{ -h - 5). 

Similarly, the jump condition is applied to momentum, and the changes in the mean 
components of velocity arc: 
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* 


AU = U{ —/i) - U{ -h - 5) and AV = V( - h ) - V{ -h - 6). 

The most distinguishing feature of the NTS model is its use of throe component 

r 

equations to define total kinetic energy (TKE) or The TKE equations, assuming 
horizontal homogeneity are: 






<0 


where pressure and water density are represented by p and p, respectively. Other terms 
appearing in these equations are those associated with the planetary rotation rate Q and 
the Coriolis parameter f defined: 

f = 2Q sin(<f>) where <p represents the latitude. 


It should be noted that when summing the contributions of these equations, all ex¬ 
pressions associated with planetar} - rotation disappear. Thus the equations can be re¬ 
duced to: 
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- 1 ¥-{(¥h¥F(¥) ; ] 


dz 


m 

(8c) 


The boundary conditions at z = 0 and at z = - /i are: 


- t/V(0) = 


Po 


u'w'(-h) = -rAU-~ 

v dt 


- v'vv'(O) < 


Po 


v'w'i -h ) = - TA v4r- 
ct 


b'w'{ 0) = + - ^ - - ' p)S b'w'( -h) = - TA5 #- 

p 0 C P Po v ' c/ 

where T is a step function defined as 1“ = 0 for ~ < 0 and T = 1 for 4r* > 0* 

c7 dt 

Garwood (1977) integrated the buoyancy equation and the momentum equations 
across the mixed layer to obtain the following bulk TKE equations. 


o{hu' 2 ) A (t » 2 . /itrv2 . „ I) 

’« 2m 3 ^ + (A U) w e + R x - — 


ct 


c(M^) 


(t v /p) 


Z) 


C7 


= 2^3 ~u, ^ (A H w e + R y~ 2 


8(hw' 2 ) 


a. 


Z) 


07 = " ^ pf + ^ /lS ( e -P)~ gh(oAT-pAS)\v t + R 2 - 


(9a) 

(9/;) 

(9c) 


New terms appearing in these equations are defined as follows. Viscous dissipation is 
D = 2tft,[(£)*'* + p 3 f/i]£ . Pressure redistribution is represented by R x , R r , and Z? 2 which 
are defined as: R x = 2m 2 {E-W)E'» ; R y = 2m 2 {E - ; R z = 2m 2 {E- W)E :2 , 

The specific heat of sea water is denoted by c P . The terms m 2 , m 2 , and p % represent 
dimensionless model constants used by Garwood to parameterize higher order processes 
into a second order closure scheme (1977). A brief description of these constants is now 
given as provided by Gallacher et al. (1983). 

The m 2 term establishes the vertical integral of the pressure redistribution as a lin¬ 
ear proportionality of the distribution of TKE. It assumes that the pres sure-strain rate 
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interaction tends to restore equal distribution of energy among the three components. 
In other words, it assumes the turbulence returns to isotropy in agreement with the 
theoretical work of Rotta (1951). The m 3 term represents the linear parameterization 
of wind shear production of TKE as a function of surface wind stress and the horizontal 
friction velocity defined in Equation 1. The combination of m, and p } shown in the 
viscous dissipation term is used to establish a parameterization of dissipation which in¬ 
corporates a Rossby number for the mixed layer.. This parameterization is another 
unique feature of the model, since it allows dissipation to be calculated in each of the 
three components and incorporates the combined effects of planetary rotation and the 
turbulent velocity relative to the length scale of the large scale turbulent flow. 

Following Garwood (1977), Equation (9) can be considered simply as the total 
vertically integrated TKE being equal to the sum of the bulk quantities of wind shear 
production, entrainment shear production, buoyant damping or production due to sur¬ 
face buoyancy, buoyant damping due to entrainment, and viscous dissipation. 

The TKE equation integrated across the entrainment zone provides the final 
equation in the model. From Eq. (4) through (6) and assuming (IV = 0) , to deepening 
of the mixed layer is 


ch _ /7 '4\'(» ') 2 E 

ci HAD + E 


( 10 ) 


The m 4 term seen in this equation is a non-dimensional model constant that represents 
the ratio of buoyancy flux to convergence of energy flux at the base of the mixed layer 
(Garwood 1977). 

Previous studies using the NTS model have assumed that entrainment shear pro¬ 
duction is negligible. Consequently, their numerical simulations omitted the entrainment 
shear production terms seen in Equations (9a) and (9b) ([(AC) 2 and (Af) 2 ]^) and the 
entrainment zone's TKE at the base of the layer shown in Equation (10). Omitting this 
process in numerical versions of the Garwood model had been an acceptable practice, 
since it facilitated computations of TKE without using the momentum equations. This 
enabled the numerical models to operate more efficiently, and it was believed by 
Garwood that any differences in the results were insignificant (Martin 19S5). 

In this investigation with the NPS model, the above processes will be included and 
the entrainment zone will be treated as having a finite thickness <5. The thickness of the 
entrainment zone is assumed to adjust so as to maintain stability wherein the product 
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of the Richardson number and the thickness of the zone (Rid) is a constant equal to the 
critical Richardson number (Garwood 1977). As shown in Figure 1, the Richardson 
number within the region -h - <5 < z < -h is treated as the critical Richardson number 
for neutral stability and is assigned a value of 1/4. The Richardson number is deter¬ 
mined by the following relationship: 


Ri- 


Otg' 


£L 

02 


Jg 4r- 

02 


A BS 


( i £)2 + ( 4 L )2 

^dz ^ dz 


(Atf + (A if 


( 11 ) 


Using the assumed value of Ri = 1/4, 6 is then obtained as a function of the changes in 
mean momentum and buoyancy across the entrainment zone: 


f(A(f+ (An 2 ](l/4) 
A B 


( 12 ) 


This completes the presentation of the theoretical description of the NPS model. 
Garwood (1976 and 1977) provides a more complete description of the theory, partic¬ 
ularly regarding the basis of many of the applicable assumptions. 
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III. 


METHOD 


A. RESPONSE TIME ANALYSIS 

This study begins by investigating the theoretical time scale for which entrainment 
shear production of TKE is significant, de Szoecke and Rhines (1976) suggested that the 
most significant contributions made by entrainment shear production occur on a time 
scale that extends from approximately 1 to 12 hours, peaking at about one half an 
inertial period at mid-latitudes after a given wind forcing event over an initially quiescent 
body of water with a linear density stratification. Additionally, the contribution of 
entrainment shear production was not shown to be a significant factor in mixed layer 
dynamics. However, their study was conducted with the Niiler model (Niiler 1975). 

The Xiilcr model differs from Garwood's model in two fundamental aspects. First, 
it establishes mean TKE as a constant function of the surface wind stress, whereas the 
NPS model computes values for mean TKE using the vertically integrated TKE 
equations. Second, it restricts dissipation to the wind-induced energy at the surface, 
rather than a more realistic estimation based on vertical and horizontal components of 
the TKE budget. Because of these differences, it was worthwhile to attempt a reprod¬ 
uction of the results obtained by de Szoeke and Rhines with the Niiler model and to 
compare them with results obtained by performing an identical experiment with the 
Garwood model. Some of the parameters in these experiments were not necessarily as¬ 
signed values identical to those used in their study, since their report did not discuss 
them. However, reasonable values approximating mid-latitude oceanic conditions should 
provide comparable results to those of de Szoeke and Rhines. 

The two models were assigned identical initial conditions and simultaneously sub¬ 
jected to wind forcing over a fluid initially at rest with a linear density stratification. It 
was found that the Garwood model initially deepened the mixed layer more slowly than 
the earlier model. The largest difference in calculated layer depth, h, occurred about 24 
hours after commencement of wind-forcing where the Garwood model computed a depth 
that was three meters less than the depth calculated by the Niiler model. As the simu¬ 
lation continued, the mixed layer depths estimated by the two models converged to 
similar solutions. These results were not surprising, because the two models differ in 
their formulation of TKE. However, it was not yet clear how these differences would 
affect the final result. 
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Forming ratios of terms used in the Niiier mode! equations, dc Szoekc ami Rhinos 
obtained results which are illustrated in Figure 2 below. 



Figure 2. Numerical solution nitli Niiier model: Results obtained by dc Szoekc 
and Rhines. (Solid curve: wind induced energy flux less dissipation; 
dashed curve: entrainment shear production; dotted curve: mean TKU.) 


The dashed curve plotted in Figure 2 represents the ratio of entrainment shear pro¬ 
duction to entrainment damping as formulated by the Niiier model. As mentioned ear¬ 
lier, entrainment shear production values are seen to provide the largest idative 
contribution to TKE between approximately 1 and 12 hours after initialization. Addi¬ 
tionally, the ratio represented by the dashed curve in the figure never exceeds a value 
of 1. 

Following their example with the Niiier model, time variant ratios of entrainment 
shear production to entrainment damping were formed for both models, and outputs 
from the simulation mentioned above were evaluated. The Niiier model computations 
indicated that the ratio of entrainment shear production to entrainment damping was 
dominant from approximately 75 minutes to 9 hours after initiation ofwind-foicing, and 
the maximum value attained by this ratio was 0.62. These results arc in gcncial agree¬ 
ment with the findings of dc Szoeke and Rhines. However, this was not the case for the 
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Garwood model. Under identical initial and forcing conditions, entrainment shear pro¬ 
duction dominated the TKE budget from 18 minutes after initialization until ab uit 12 
hours later. Additionally, the maximum value attained by the ratio was 1.52, which was 
over twice the highest value computed by the Niiler model. 

The question arises as to why these results differ. The results are inheicntly i elated 
to the different deepening rates of the models. The latio used Tor computing these values 
from the Niiler model was: 

entrainment shear production 2/ _2 m(l— cos(/7)) 
cntiammcnt damping -L/V 7 /, 11 

whcic /i, rcpicscnls the mixed layer depth computed by the Niiler model. A complete 
dcsciiption of the individual terms in the above expression is provided by Niilci (1975). 
Similarly, the Garwood model's computation for the same ratio can be cxpicsscd as: 

2/~V(l- cos {ft)) 

y /V’/iJ-Ev® 

where /i, is the layer depth calculated by the Garwood model. 

For this particular simulation 6'w'(0) equaled zero, so the two models performed 
• ' identical computation of this ratio. However, since the models yielded different val- 
"cs for mixed layer depth, their outputs differed by a factor of (/fa/A,) 4 . Ihc magnitude 
of this factor happened to be greatest duiing the first inertial period when shear pro¬ 
duction played its most active role. As an illustration, consider the following values ic- 
tuincd by the models two hours after initiating wind forcing. 

Niilcr's ratio = 0.57 Garwood’s latio *= 1.25 
/i, = 1148.2 cm /ij = 943.3 cm 

The iatio of the fust two values (0.57/1.25) has a value of 0.456, which is equal to the 
value of (hjIliiY • 

Ihc appaicnt significance of entrainment shear production is thus shown to vary 
between these two models according to their differences in computing the ’1 KE and in 
accounting foi dissipation. Since the Garwood model uses a moic compiehcnsive 
paiametciization, these icsults suggest that entrainment shear production plays a moie 
significant tolc in the ’1 K12 budget than may have been predicted by cailici mixed layer 
models. 
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Although de Szoeke and Rhines (1976) suggested that this process is associated 
with a particular time scale, it is noted that their study was restricted to a column of 
water having a linear density stratification initially. A more accurate assessment of the 
significance of entrainment shear production migtit be obtained by evaluating the proc¬ 
ess with respect to mixed layer depth. The largest relative difference between the layer 
depths predicted by the Garwood and Niiler models (A, - h 2 ) occurred when the layer 
was shallow, thus (/i 2 //i,) 4 attained a high value. Additionally, the value oi h 4 appears in 
the denominator of the expressions used in our analysis, which creates an inverse re¬ 
lationship between the mixed layer depth and the contribution of entrainment shear 
production. These results suggest that the process of entrainment shear production is 
most significant when the mixed layer is shallow, implying that simulations of the upper 
ocean at Ocean Station Papa should be most affected by this mixing process in the 
summer and early fall seasons. This investigation now continues with a series of tests 
of the NPS model. 

B. PRELIMINARY ADJUSTMENTS 

The version of the NPS model including entrainment shear production used in this 
study is in the developmental stage, and a few changes and corrections in its numerical 
code were necessary before proper testing could occur. Each coding change required 
verification, and meant that all previous results were suspect. It is believed, however, 
that the NPS model with entrainment shear production added received adequate testing 
during the course of this investigation to render it suitable for future studies of air-sea 
interaction or for embedding into general circulation models such as was done with an 
earlier version by Adamec et al. (19S1). 

The XPS model with entrainment shear production is identical to the version 
without entrainment shear production, except that it performs additional computations 
to simulate this additional source of TKE within a zone of finite thickness below the well 
mixed surface layer.. Both versions were verified using Ocean Station Papa data ex¬ 
tending from January 1,1961, to December 26, 1969. Each run was initialized with ar¬ 
chived BT and bucket SST data. In each simulation, the temperature profile recorded 
nearest in time to the simulated initiation time provided the initial conditions of the 
thermal structure of the upper 200 meters of the ocean at Station Papa. To perform 
operations with a vertical resolution of I m, the BT data were linearly interpolated be¬ 
tween the temperature recordings at five meter intervals. A time step of 1 hr was used, 
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and atmospheric forcing was computed and interpolated using thrcc-hourly 
meteorological observations and bulk flux formulas. 

Model predictions were available at one hour intervals, but for investigating long 
term simulations (seasonal and annual) it seemed adequate to use data output every 
third hour. Model SST prediction was chosen as the primary variable for comparing 
model performance against observed bucket SST. This had the advantage of simplifying 
comparisons of large sets of data. Some difficulties arose because of this choice, and 
these will be discussed later. 

Having interpolated the initial temperature profile to a resolution of 1 m, the NFS 
model then computed successive temperature values for each grid point extending from 
the surface to a depth of 200 m. For initialization purposes, the mixed layer depth was 
defined as the shallowest depth at which the temperature was 0.2 C less than the SST. 
The details of the entrainment zone numerical scheme used in the NFS model are pro¬ 
vided by Adcmcc ct al. (1981). 

Since the original numerical version of the model had been modified to include 
simulation of entrainment shear production, it was appropriate to verify the conserva¬ 
tion of heat. For this purpose, 


-1 200 


heat change* j-f QJz)<!z Qj{z)dz oc]T( 7) - T f )j 

J-h-6 J j -1 


where the time variant heat content ( Q ) at each grid point in the water column is pro¬ 
portional to the change in temperature at the appropriate level (J), and initial and final 
values are denoted by subscripts i and /, respectively. Choosing reasonable values for 
model parameters, several 96 hour simulations were conducted with the new version of 
the model. For each simulation, the above computation was performed, and typically 
returned error values of less than 10' ,: C. These small values could be attributed to ma¬ 
chine round-olf errors, so it seemed reasonable to assume that the model was operating 
within the constraints of conservation of heat and potential energy, 

The earlier discussion regarding the theoretical time response of entrainment shear 
production suggests that the contribution of this process should vary depending on the 
depth of the mixed layer. This led to an examination of the seasonal performance of the 
NFS model, 

The version of the model without this additional process was integrated for 90 day 
simulations corresponding to seasons from January to March, April to June, and so on 
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for the period extending from January 1,1960, to December 30, 1968. These simulations 
yielded three-hourly outputs which included model predictions ( T moi ) and bucket obser¬ 
vations { T rf,) of SST. From these data seasonal values of root-mean-squarc (rms) error 
were examined where this error was defined as: 


1 V 2 

rms = “ / (<57mod ” 6 T ° b ^ 

/=! 

New terms appearing in this expression are A', the number of three-hourly data points, 
and ST mai and 5T oks defined as the difference between instantaneous and mean values of 
model and observed SST, respectively. It appeared that the performance of this version 
of the model was best in the January to March time-frame and poorest in the October 
to December period. Typical rms error values for the winter season were appoximatelv 
0.2 C 2 . while those for the Fall season were 0.7 O with a variability of about + 0.16 O 
for each. 

Before concluding that these results were significant and that model performance 
could be categorized by season, it was noted that the losults were influenced significantly 
by the particular day selected for model initialization. For example, choosing to start a 
90 day simulation on September 28. 1966, versus September 29, 1966, yielded error re¬ 
sults that were significantly different as illustrated below: 

• From 28 September: mean error = 1.105 C and rms error = 0.55S O 

• From 29 September: mean error = 0.933 C and rms error = 0.435 O 

where mean error was defined as the average of the difference (7' mod — T ci>! ) of the three- 
hourly outputs for the simulation. Thus the "seasonal" performance of the model could 
possibly be attributed to this sensitivity to initialization time. Earlier research by 
Warrenfeltz (1980) examined the sensitivity of this model to changes in intitial condi¬ 
tions and noted, in particular, significant differences in results for 15 day simulations 
initiated with different conditions during the fall deepening of the mixed layer. Thus it 
seemed that this approach to evaluating the seasonal performance of the model would 
not yield statistically significant results. 

At this point a different method of studying the long term effects of entrainment 
shear production was adopted. The length of model simulations was extended to 360 
days, and model initialization was fixed at 1 January for each year. Initializing the 
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model in the winter was preferred, since the mixed layer was deep and not in transition. 
Additionally, it facilitated making general observations of inter-annual variability and 
comparing model results with previous studies such as Martin (1985) and Gaspar(1988). 
The seasonal characteristics of model performance could be examined for tendencies to 
degrade or improve with time based on cumulative errors brought forward from a pre¬ 
vious season. 

C. TUNING THE MODEL 

Before obtaining meaningful results from simulations with the NPS model, some 
criteria needed to be established for selecting values of the non-dimensional model con¬ 
stants (eg., nt 3 , /> 3 ) mentioned in Equation 9 and its subsequent discussion. Although 
these constants should theoretically be universal in their applicability, the assumptions 
used in developing this second order closure one-dimensional bulk model have neces¬ 
sarily neglected a complete description of the physical processes involved in the dynam¬ 
ics of the upper ocean. As shown by Gallacher et al. (1983), the values assigned to sonic 
of these constants have varied considerably in previous studies. 

Since this investigation would mainly use SST errors as the criteria for comparing 
model results, it was decided to tune the models by seeking an optimal combination of 
model constants for reducing SST error. The NTS model without entrainment shear 
production had already been tuned to optimize its performance against data at Ocean 
Station Papa for certain model simulations, but not specifically for SST prediction cor¬ 
responding to a given year. Most model constants were assigned values that had been 
found acceptable in previous studies. These included the following: tn { — 1 , m 2 — 1 , and 
m t = 1. The model constants which seemed to have the greatest influence on model re¬ 
sults and which would have been most affected by including the additional process of 
entrainment shear production were /n 3 and p y For these reasons, the following approach 
was taken to tune the model for this study. 

First, a four year period from 1966 to 1969 was selected to represent a series of 
w'ell-observed years at Ocean Station Papa, and an extensive scries of 360 day simu¬ 
lations was conducted for each year with each version of the model. Every run w r as ini¬ 
tialized with the nearest archived temperature profile corresponding to January 1 for 
each year. The only differences between each simulation were the values selected for 
nt 3 and p y 

Next a 2-D array of rms error was formed for each of the two model versions as 
a function of the independent variables, m 3 and p } As these arrays were computed and 
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plotted, a minimum rms error was sought for each year along with its cot responding 
values ofmj and p y Although some inter-annual variability existed in the best values for 
these two constants, the overall pattern of the array fields of rms error constructed for 
a given version of the model maintained a consistent pattern in the four years investi¬ 
gated. Tor the version without entrainment the optimal values for m 3 and p 3 were found 
to be in the vicinity of two to four. The results with the version with entrainment shear 
production suggested that p 3 should have a value of approximately one to three and m 3 
should be about equal to one. 

Having found that the best values of the model constants varied from year to year, 
it was decided to tunc the model to achieve optimal results for a single rcpicscntativc 
year. The year 1966 was chosen as the base year, since the array fields generated by both 
versions of the model for that year appeared closest to representing the results obtained 
for the four years examined Using data from 1966 each version of the model was re¬ 
peatedly integrated over the year to provide optimum values of the constants to an ac¬ 
curacy of the nearest one hundredth. This approach led to the following values for the 
base year: 

• Model version without entrainment shear production: p 3 - 2.56 and m 3 = 2.98 

• Model version with entrainment shear production: p 3 «* 1.44 and m 3 = 0.90 

It is noted that since the two versions of the model are each attempting to represent 
higher order processes with a different set of physical processes, the optimal values for 
the model constants should differ. Additionally, the version with entrainment shear 
production is more sensitive to the values selected for these model constants. This sen¬ 
sitivity is illustrated by the array fields seen in Figure 3 on page 20. The greater sensi¬ 
tivity of the newer version of the model may suggest that it is more realistically 
representing the dynamics of the upper ocean. 

Inter-annual variability causes the minima shown in Figure 3 to shift, particularly 
along the axis representing values for m y Observations taken from 1966 through 1969 
show that the minima migrate toward lower values of m 3 when Ocean Station Papa da'a 
contains SST values which are higher than those for 1966. As shown by the solid curve 
in Fig. 2, wind shear provides the dominant contribution to the TKE budget, so attempts 
to match observed SST values with the model inevitably involved adjusting the 
parameterization of the wind forcing term. This was accomplished by changing the 
magnitude of m y The combination of this modeled feature (which was observed in both 
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versions) and the increased sensitivity of the newer version to values of m 3 adversely af¬ 
fects efforts to find universal values for the model constants. However, for the purposes 
of this study, the inter-annual variability is not an obstacle to examining the significance 
of entrainment shear production. 


19 




Model version without en 


o 

t o."' 

L»J 



Model version with entrai 

Figure 3. Error arrays: Plots of rms ' 
dieted SST using data from tin 














IV. RESULTS 


A. OVERVIEW 

Since the model had been tuned to optimize annual simulations of SST at Ocean 
Station Papa in 1966, this presentation of results will highlight SST simulations during 
that particular year. Model SST performance against data from other years at the sta¬ 
tion will be discussed briefly, as will simulations of the depth of the mixed layer. These 
will be followed by introducing the improved vertical representation of the upper oceanic 
temperature profile that is provided by adding the treatment of entrainment shear pro¬ 
duction in a finite zone. Finally, results will be presented on the significance of 
entrainment shear production in mid-year simulations of the mixed layer in mid- 
latitudes. 

B. SST SIMULATIONS 

The earlier studies by Martin (1985) and Gaspar (1988) showed that annual simu¬ 
lations of SST bv bulk mixed laver models tend to demonstrate two intrinsic features. 

r •> 

The first of these features was pointed out by Martin who noted that these models tend 
to yield less variability in simulations of SST than is observed in the data on the synoptic 
time scale during the fall and winter. In addition, he noted that these models yield more 
variability than is observed in the data on the synoptic time scale during the spring and 
summer. The second feature was discussed by Gaspar who noted that these models tend 
to over-estimate SST values (while under-predicting h ) in the summer and to under¬ 
estimate SST values (while over-predicting h) in the autumn. Gaspar observed that ad¬ 
justing the model parameters to correct the SST errors for one season increased the 
errors in the other.. 

Having tuned the NPS model to reduce annual rms errors in the estimates of SST 
during 1966, the "optimal'' results for that particular year were first examined in light of 
the features noted in the studies by Gaspar and Martin. On the annual time scale results 
from the two versions appeared almost identical, and the principal features in the fol¬ 
lowing discussion apply to both versions except where noted otherwise. Hereafter, the 
version of the model that includes the additional process of entrainment will be referred 
to as the "enhanced version". 

The results of the 1966 simulation of SST by the enhanced version are presented in 
Figure 4 on page 25. The behavior noted earlier by Martin is evident, although the pe- 
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riods studied here are not exactly coincident with those of Martin. The tendency of the 
model to demonstrate less svnoptic-scale variability in SST values than is recorded in the 
observations is shown in Figure 4 to occur from January till late April and then again 
after late September. The opposite tendency is seen in the warmer months from late 
April until late September. In Figure 4 two prominent examples of these two tendencies 
are shown by the different behaviors of the two curves on January 22 and on June 21. 

The tendency of the model to exhibit less synoptic-scale variability in the cooler 
months appears to be associated with the relatively larger depths of the principal (or 
seasonal) thermocline. With deep mixed layers it is difficult, using only heat and mo¬ 
mentum fluxes, to adequately simulate synoptic-scale variability in an inhomogeneous 
upper ocean. Perhaps other factors (such as advcction or upwelling) that contribute to 
the heat content of the upper ocean provide an increased relative contribution when the 
mixed layer depth is relatively larger. Since the model has been tuned to best represent 
the annual evolution of SST, it attempts to compensate for these missing physical 
processes by over emphasizing the roles played by the fluxes of heat and momentum. 
This over-emphasis would explain the model behavior observed in the warmer months 
when the principal thermocline has shoaled, and the uniform thermal characteristics of 
the mixed layer are more easily influenced by surface heating and wind-mixing. This 
intrinsic feature of the XPS model was observed in simulations against data from other 
years at Ocean Station Papa, and it could not be removed by adjusting the model con¬ 
stants. 

Gaspar's observation that bulk mixed layer models over-estimate SST in the summer 
and under-estimate its value in the autumn was not demonstrated by the NPS model 
with the constants tuned optimally for 1966. As Figure 4 shows, the model slightly 
under-predicted the observed SST throughout the year, except for about a three week 
period in the late autumn. To further examine the general tendency of the model to 
under-estimate SST in an annual simulation, additional runs were performed with two 
significant variations. First, to determine if this tendency was unique to 1966, the model 
was tested aaainst data from nine other years usine the same model constants. Addi- 
tionally, to see if the tendency was a by-product of tuning the model constants to reduce 
rms errors in SST values, tests were made against 1966 data with significant variations 
in the values assigned to m ; and p y 

In testing the NPS model against data from other years, an error analysis was con¬ 
ducted for each 360 day simulation. These errors were compiled to provide a general 
evaluation of the model's performance against data extending from 1960 to 1969. Errors 
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for the annual simulations are provided in Table 1, and it is seen that with the optimal 
model constants for 1966, the model always tended to under-estimate SST. Also indi¬ 
cated in the table is a significant amount of inter-annual variability, with rms errors 
ranging from 0.29 O to 1.24 O and mean errors recorded between -0.12 C and -1.7 C. 
To better understand the diversity shown in these results, the last five years of this period 
were examined in detail and will provide the basis used in the following general de¬ 
scription of annual simulations of SST with the NPS model. 

From January to April, SST prediction tended to be quite accurate, with errors sel¬ 
dom outside of +0.5 C and -0.5 C. From May through October, the model always 
tended to under-estimate SST with errors varying greatly from year to year. The under¬ 
prediction of SST in the late summer was typically -1.0 C, but for 1965 these errors fre¬ 
quently surpassed -3.0 C. The model tendency during this period was the principal 
reason for the mean errors recorded in Table 1, and especially for the unacceptably large 
value shown for 1965. The last two months of the year returned errors ranging from 
about -1.6 C to + 1.4 C, and no significant trend was noted except that the performance 
near the end of the year was related to how well the model had predicted SST values 
through the late summer period. Thus, the tendency of the model to under-estimate SST 
values as shown in Figure 4 is correct only in that this tendency is confined to the 
warmest months of the year. Also the magnitude of error illustrated in Figure 4 is not 
necessarily representative of other years. 


Table 1. ANNUAL SST ERRORS. 
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The next series of runs conducted with the model were performed to examine how 
tuning the model constants might affect the model's tendency to undcr-cstimate SST 
throughout most of 1966. Especially desirable was evidence that the tendency could be 
reduced in the warmest months without adversely affecting model performance in the 
other times of the year. It was found that either increasing the value of p 3 or decreasing 
the value of w 3 eliminated the tendency of the model to under-estimate SST throughout 
the year. Increasing the value of/> 3 represented an increase in the dissipation rate of TKE 
in the mixed layer, and reducing the value of m 3 represented a decrease in the generation 
of TKE by wind-forcing. Both of these changes resulted in less mean TKE in the mixed 
layer, permitting the SST predictions to reach higher values. Although for future studies 
some combination of changes to both of these model constants would be best, it was 
noted that of the two, increasing the value of /> 3 had less of an adverse effect on simu¬ 
lations of the cooler periods of the year. A rough estimate from these runs is that it 
should be assigned a 50 per cent larger value than that assigned to each version of the 
model in this study. 

A noteworthy feature shown in Figure 4 is the abnormally large SST observation 
recorded on January 22. The value shown is about 7.7 C and was obtained by bucket 
measurement. This value did not appear in the mechanical BT data which consistently 
showed near surface temperatures of less than 6.0 C from January 21 to January 23. 
Although the larger value recorded in the bucket data is considered spurious, no attempt 
was made to remove it or any other data from either the simulations or the error ana¬ 
lyses. 

C. UPPER OCEAN THERMAL STRUCTURE SIMULATION 

This aspect of model results will highlight some of the general tendencies exhibited 
in upper ocean simulations with the two versions of the NPS model. As was mentioned 
in Chapter 3, model initialization could play a significant role in the outcome of a given 
simulation. Initialization with BT profiles from either of two successive dates could 
greatly alter the results. The adverse factors contributing to an erroneous BT are nu¬ 
merous, including navigational error of platforms, calibration of BT's, precision of tem¬ 
perature recordings (~ 0.1C), and a vertical resolution of only five meters. These same 
factors contribute adversely to model verification on any given date throughout a simu¬ 
lation. Although these factors make it difficult to assess errors in the model, they have 
the same effect on results from both versions. A brief description of tendencies common 
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Figure 4. Model Estimated SST: Values yielded from 360 day simulation for 
1966 (Solid line: model SST, dashed line: observed SST). 
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to both versions will be presented, focusing on some of these problems. This will be 
followed by a presentation of differences exhibited between the two versions of the 
model, highlighting their skill during the warmer months and their physical represen¬ 
tation of the upper ocean. 

The ability of the NTS model to simulate the annual evolution of the thermocline 
is presented in Figures 5 and 6. The asterisks in each figure depict data points from BT 
data archived at five meter intervals, while the solid and dotted lines represent the tem¬ 
perature profiles predicted by the enhanced and earlier versions of the model, respec¬ 
tively. As illustrated in Figure 5 on page 28, a problem occurred in the initialization of 
the model. The 21 January profile shows an observed mixed layer isotherm of 5.8 C. 
Both versions showed a -0.4 C error in the temperature of the layer three weeks after 
initialization. This represents a very large amount of heat, since the layer depth was 
approximately 135 meters. Two days later, the 5.S C isotherm extended only to a depth 
of 100 meters: the next 35 meters had a lower value of 5.7 C. On 30 January, one week 
later, the 5.9 C isotherm was found to exist at these depths, and to further complicate 
matters, the temperature at the surface had decreased to 5.7 C. The SST for 21 and 23 
January was 5.8 C, which was the value of the 30 January isotherm between 15 and 30 
meters. This evolution of temperature structure in the upper ocean cannot be explained 
by one-dimensional processes. Physically speaking, this sequence of BT profiles would 
imply heating from the lower ocean, which is not possible. This phenomenon is best 
explained by one of three possibilities: ship movement, advection, or observational er¬ 
ror. In other words, either these BT data were collected from different water masses, 
or temperature variations of ± 0.1C should be considered insignificant. 

The problems associated with the temperature profiles observed in late January had 
disappeared by late March, as shown by the 26 March profiles in Figure 5. It is noted 
that the estimated depth of the mixed layer on 26 March is the depth at which the model 
estimate of temperature is approximately 0.6 C (for the earlier version) to 0.7 C (for the 
enhanced version) less than the simulated SST. Also the model temperature profiles, 
including the depth of the mixed layer, are in good agreement with the observed BT data 
at that time of year. This is fortuitous, since it lends more credibility to a discussion of 
model performance during the warm summer months than would have been possible 
after the initialization problems observed in January. 

As seen in Figures 5 and 6, the general tendency of the model was to under-estimate 
the depth of the mixed layer throughout the summer months and into the fall. In view 
of the tendencies observed in the previous discussion of model prediction of SST, this 


26 




tendency to under-estimate the depth of the layer is significant for two reasons. First, 
since the model tended to under-estimate SST throughout the warr ter months, this ad¬ 
ditional tendency implies either an inaccurate heat llux at the surface or an inaccurate 
amount of mixing in the upper ocean. This tendency was also demonstrated by the 
one-diinensional models examined by Gaspar (1988), and further research will be re¬ 
quired to select the correct reason. The other significant aspect of the model's tendency 
to under-estimate layer depth during the warmer months is the synoptic-scale variability 
problem noted in the SST simulations. A smaller layer, as predicted by the model in the 
warmer months, would be more responsive to synoptic-scale forcing. Since the model 
exhibited more synoptic-scale variability in the summer months than was observed for 
SST values, the parameterization of momentum and heat forcing may not be as poor 
as was suggested by the results in those simulations. The thinner layers predicted by the 
model should be more responsive to synoptic-scale forcing than the deeper layers indi¬ 
cated by observations. Thus the problems observed for SST estimations arc probably 
related to a long-term error or bias in the net heat flux imposed at the surface. This is 
not a model problem but a boundary 7 condition problem. 

An interesting feature observed in the simulated upper ocean thermal structure was 
the evolution of a thermocline at about 60 meters depth. As Figure 6 shows, this feature 
represented a mean temperature change of approximately 1.2 C, and corresponded to 
significant mixing in the upper ocean sometime after 26 March. On that date, as seen 
in Figure 5, the SST was 5.9 C and the mean temperature of the layer was 5.8 C. 
Meteorological data shows that a large storm passed by Station Papa between 12 and 
15 May, and as Figure 4 shows, the model estimated SST fer those dates was about 6.1 
C on May 12 and 5.8 C on May 15. Thus the thermocline feature was created by the 
model in response to wind forcing during that storm, and it remained in the temperature 
profile until its removal by mixing during the fall deepening. This identical feature did 
not exist in observed BT profiles examined in this study. Perhaps the gradient in the 
observed data at about 70 meters depth in the 27 May profile corresponds to the same 
storm. It is noted that the observed SST in the isotherm at the base of the gradient is 
6.3 C which corresponds to the observed SST value on 15 May. The observed gradient 
at about 70 meters persisted with some variability (perhaps due to instrument precision 
or ship movement). However, through the summer months, unlike the model 
thermocline at 55 meters, the observed gradient exhibited much less intensity (less 
A7‘/Az) by September. This suggests that diffusion is greater in the upper ocean than it 
is in the model. Further research is needed to understand how to improve 
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Figure 5. Temperature profiles: Solid line: enhanced version estimate, dotted 
line: earlier version estimate, asterisks: observed BT data. 
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Figure 6. Temperature profiles: Solid line: enhanced version estimate, dotted 
line: earlier version estimate, asterisks: observed BT data. 
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parameterization of heat conduction and ambient diffusion below the surface mixed 
layer and entrainment zone. 

D. EXAMINING ENTRAINMENT SHEAR PRODUCTION 

The annual simulations of SST by the two versions of the NPS model did not differ 
significantly, since each version had been tuned to optimally reproduce the observed SST 
values. However, there were notable differences in the temperature profiles yielded by 
the two versions as shown in Figures 5 and 6. Although these differences might be at¬ 
tributed to entrainment shear production, it must be pointed out that to some extent, 
differences had been introduced by using different values of p 3 and m 3 for the two ver¬ 
sions. For these reasons, two approaches were adopted to evaluate the significance of 
entrainment shear production. The first approach was to observe temperature profiles 
yielded by the two versions after they had been modified to remove the biases created 
by independently tuning them to reduce SST errors. This was easily accomplished by 
assigning a value of 2.0 to both m 3 and p 3 for the two versions. The second approach 
involved computing a time series of the differences between SST values yielded by the 
two versions and examining these differences for a relationship to features observed in 
a time series of the atmospheric forcing. 

The response analysis conducted earlier suggested that entrainment shear pro¬ 
duction should be most significant during the warmer months when the mixed layer is 
shallow. It will now be shown that the results obtained by equalizing the model con¬ 
stants for the two versions suggest that this tendency exists, but perhaps is only mani¬ 
fested in features very near the surface, especially in estimates of SST. 

The temperature profiles hindcastcd by the two versions of the model on 21 January 
were nearly identical. However, as shown in P r igure 7 on page 32, the two versions be¬ 
gan to yield slightly different results by late March. By late July the results were 
measureably different. The summer profile shows that significant differences occurred 
as a conseqence of each development of a major thermocline. The enhanced version 
created the mid-May-storm thermocline at a depth of about 68 meters, whereas the 
version without entrainment shear production established the same thermocline at about 
52 meters. Furthermore, the temperature at the top of this steep gradient was about 3.5 
C cooler in the enhanced version's results, corresponding to the increased mixing which 
had mixed the layer to a greater depth. This phenomenon was evidently repeated for 
each major storm period at the station, as evidenced by the differences in the two other 
prominent thermoclines shown for 24 July. One major thermocline is shown at about 


30 




32 meters in the enhanced version, contrasted with a depth of about 22 meters in the 
earlier version. The other major thermocline is shown near the surface for each version. 
It is interesting to note that between the two versions, the smaller the difference in the 
depths shown for a particular thermocline, the larger the difference in the total temper¬ 
ature change (A7). Thus the two versions returned only a 0.4 C difference in SST after 
the mid-Mav storm, but show a 3.6 C difference on 24 July. These results are consistent, 
since they imply that the treatment of the total heat budget in the upper ocean by the 
two versions is quantitatively the same. The differences could thus be attributed to the 
increased mixing in the enhanced model provided by the additional process of 
entrainment shear production. However, since differences in the results began to appear 
as early as March, it remains to be demonstrated that the process is more significant 
during the summer months. To this point, the summer months show only a more pro¬ 
nounced difference in SST between the two versions. 

The second approach to examining the significance of entrainment shear production 
involved comparing the results of the two optimally tuned versions in view of atmo¬ 
spheric forcing. Two aspects of atmospheric forcing were studied to gain an under¬ 
standing of the roles of the fluxes at the ocean surface in entrainment shear production. 
These aspects were wind stress (r) and net downward heat flux (Q 0 ), plots of which are 
provided in Figures 8 and 9, respectively. These quantities were calculated using stand¬ 
ard bulk aerodynamic formulae with a non-dimensional wind drag coefficient (C„) of 
1.3 x 10 3 . The employment of these formulae at Station Papa was well described by 
Raney (1977), and the values assigned to constants used in the formulae were adopted 
from that study. To compare the results of the two model versions, three-hourly outputs 
of SST were used to compute values of the difference in SST provided by the two model 
versions. A time series plot of these differences is provided in Figure 10 where each 
three-hourly data point represents a particular SST value returned by the enhanced ver¬ 
sion minus the corresponding value computed by the earlier version. The plot in Figure 
10 has been partitioned into eight periods to facilitate a discussion relative to informa¬ 
tion provided by the flux plots shown in Figures 8 and 9. 

The first period (/) shown in the SST difference plot extended from initialization 
until a rapid succession of moderate wind-forcing events in the early summer beginning 
on about day 172. During this period the SST differences between the two versions were 
insignificant, although some variations began to apperr associated with spring storms. 
The second period {IT) was characterized by the establishment of a highly variable dif¬ 
ference between the two versions between days 172 and 232. During period II the mean 
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difference between the two versions was about -0.25 C . On about day 232 a significant 
shift occurred in the mean difference in SST values, marking the onset of period III. This 
period lasted until day 259 and was characterized by warmer SST results from the en¬ 
hanced version with a mean difference of about + 0.08 C . Period IV began with a sharp 
increase in the SST difference to almost +0.55 C on day 260, followed by a nearly 
monotonic decrease to about +0.15 C on day 269. In the following period (l 7 ) this dif¬ 
ference was maintained with very little variability. Period 17 commenced on about day 
280 and showed another almost monotonic decrease in the SST difference that extended 
to about day 290. This was followed by two other periods (VII and VIII) that each 
showed a gradual increase in the SST difference from about -0.2 C to about -0.3 C. 
These last two periods were separated by a significant event on day 332 that caused the 
difference to drop by just over 0.1 C. 

Before relating the SST differences shown in Figure 10 to the atmospheric forcing 
data, a few general features shown in the figure should be discussed along with a review 
of the conditions applicable to a hypothesis regarding the significance of entrainment 
shear production. In general, it is noted that the variability in SST differences was most 
prominent during period III and least significant during the first half of period I and 
throughout periods 17/ and VIII. These periods correspond to a deep principal 
thermocline and show that if entrainment shear production was making a contribution 
to the TKE budget, it was not manifested in a change in SST. However, low variability 
in SST differences was also exhibited during period V, a time when there was a shallow 
principal thermocline. The combination of these features suggests that the existence of 
a shallow principal thermocline is a necessary, but not sufficient, condition for 
entrainment shear production to become significant, at least near the sea surface,, Ad¬ 
ditionally, the results from the investigation of the upper oceanic thermal structure sug¬ 
gested that differences due to entrainment shear production occurred as a consequence 
of significant mixing events, as was illustrated in Figure 7. Furthermore, the results from 
the response analysis conducted earlier suggested that the process should affect the up¬ 
per ocean most, not only when these conditions are met, but also only when the ocean 
is at rest immediately prior to a wind-forcing event. In other words, the results estab¬ 
lished thus far suggest that entrainment shear production should be most significant 
under the following conditions: 

• A relatively shallow principal thermocline. 

• A shallow mixed layer, or none at all. 
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• Initially, a period of weak or absent wind-forcing. 

• A sudden, rapid increase in wind-forcing sufficient enough to deepen or create a 
mixed-layer. 

To physically isolate these conditions at Station Papa, the data were searched for 
times when the following features were exhibited: 

• A shallow principal thermocline. 

• A recent period of large downward net heat fluxes. 

• A recent period of weak wind-forcing. 

• A sudden, rapid increase in wind-forcing, particularly one that is not significantly 
prolonged. 

The persistence restriction ("not significantly prolonged") was added to the wind-forcing 
condition as an attempt to isolate a period from conditions created by a difference in the 
parameterization of wind shear production between the two model versions. The value 
assigned to the m 3 constant in the earlier version was over three times larger (2.98 versus 
0.90). Thus sustained wind conditions would prevent isolation of the effect of 
entrainment shear production. 

Somewhat arbitrarily, these conditions identified for investigating entrainment shear 
production were assigned threshold criteria. By no means were these criteria considered 
universal. They were merely introduced to provide a convenient means of analyzing the 
data., The first four of these threshold criteria were used to examine the time scries plots 
in Figures 8 and 9 and are described as follows. 

The first criterion, shallow principal thermocline, was defined to be at a depth above 
-50 meters, as estimated by the enhanced model. The next criterion, a period of large 
downward heat flux, was defined as cither a 24 hour period during which the peak net 
downward heat flux surpassed 400 fl'/m 2 or a 48 hour period during which it reached 350 
Wjm : for two consecutive diurnal cycles. Since z has been defined as positive upward, 
these values are negative in Figure 9. The third criterion, a a period of weak wind¬ 
forcing, was defined as a period of 24 hours or longer having r values less than 0.02 
dynes I cm 2 . This corresponds to a wind speed of less than I m/s at a height of 10 meters. 
The fourth criterion, a wind-forcing event likely to provide significant mixing by 
entrainment shear production, but isolated from sustained wind-forced events was de¬ 
fined as one in which t values rapidly rose to greater than 1.7 dynes!an 2 , but were not 
sustained above that value for greater than 24 hours. This criterion was designed to 


34 







provide sufficient momentum for entrainment shear production, while avoiding the bias 
created by the larger m 3 term used in the earlier model. 

Having established the previous four criteria for investigating the time series shown 
* in Figures 8 and 9 for conditions likely to provide significant mixing by entrainment 

shear production, a final threshold criterion was established. This fifth criterion was 
, used to determine if a significant contribution to the TKE budget had been made by 

entrainment shear production. It was defined as a greater than .2 C negative SST change 
in the difference in SST estimates between the two model versions occurring within a 
time period of 18 hours (corresponding to about one inertial period). To illustrate this 
definition, the reader’s attention is invited to Figure 10. In this figure, the difference in 
SST values estimated by the two versions made five negative deviations that surpassed 
this threshold criterion. These each occurred within two day periods corresponding to 
days 172 to 173, 186 to 187, 194 to 195, 204 to 205, and 232 to 233. Other less prominent 
negative deviations are shown in this time series, and each of these were examined in 
view of the principal thermocline depth and forcing criteria discussed earlier. 

The results of this examination are presented in Table 2, where the five threshold 

criteria are represented by SHPT (shallow principal thermocline), Q-l (recent warming), 

- TAU-1 (initially calm winds), TAU-2 (proper wind-forcing), and SESP (significant 

/— 

entrainment shear production), respectively. The appearance of a check mark ( v ) in 
any of the criteria columns in Table 2 indicates that the threshold criteria were met, and 
the final column in Table 2 shows verification of the hypothesis. It should be noted that 
the criteria established for this hypothesis are ad hoc. However, the qualitative results 
illustrated in Table 2 tend to verify the hypothesis. 

The only instance found where the hypothesis was not verified in accordance with 
the criteria was during a period in early July, days 186 to 187 in Table 2. Although the 
peak net downward heat flux was never very large during the two days prior to this pe¬ 
riod, it reached a very large value on day 187, and a thin layer associated with the diurnal 
cycle may have formed prior to increased wind-forcing. Also, it is noted from Figure 8 
that the winds had been weak for several days immediately prior to this event. These 
conditions also would have been favorable for a shallow mixed layer to have existed on 
day 187. 

These results suggest that entrainment shear production can indeed play a signif¬ 
icant role in the TKE budget, but only when the appropriate conditions have been sat¬ 
isfied. Furthermore, these conditions are not limited to hypothetical numerical 
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experiments, but may frequently occur in any region of variable wind forcing and shal¬ 
low mixed layers. 
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Table 2. EFFECT OF ENTRAINMENT SHEAR PRODUCTION: Symbols and 


abbreviations are explained in the text. 
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Two other sets of interesting features appear in the plot of SST differences shown 
in Figure 10. These are the significant positive and negative shifts in the SST differences 
occurring at various times throughout the year. 1 he positive shifts occurred near the 
end of period 11 on day 226. in the beginning of period IV on day 259, and between pc- 
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riods VII and VIII on day 331. Examination of the wind-forcing data in Figure 8 shows 
that each of these positive shifts was associated with either a very strong wind-forcing 
event or persistent, strong wind forcing. Either of these conditions caused significant 
wind shear production. This tended to bias the mixing differences between the two model 
versions, with larger TKE (and lower SST) estimates returned by the earlier version due 
to the higher value assigned to its wind-forcing parameter, m } . Also, these events were 
not immediately preceded by a sufficient lull to permit entrainment shear production to 
significantly alter the SST. It should be noted, however, that these conditions had less 
effect when the layer deepened, as illustated by the storm on day 291. Very intense 
forcing was required to accomplish the shift on day 332. 

The significant negative shifts are shown in Figure 10 to have occurred on three 
occasions in 1966. The first occurred during period II between days 172 and 177. The 
second and third shifts occurred during periods IV and VI. Each of these periods in¬ 
cludes a succession of small entrainment shear production events. These results suggest 
that the contribution of entrainment shear production to the TKE budget is not only 
determined by the intensity of a given entrainment shear production event, but also 
upon the number of events. 

Upon examining the comparisons in this study (one with optimally tuned versions 
of the model and one with identical model constants for the two versions), one final re¬ 
sult was obtained. This was an improved vertical representation of the physical prop¬ 
erties in the upper ocean by inclusion of entrainment shear production. A typical 
temperature profile is provided in Figure 11 showing the computed and observed tem¬ 
perature profiles for May 27< As seen in the figure, the enhanced version shows a finite 
entrainment zone at the base of the layer that reduces the temperature gradient in com¬ 
parison with that predicted by the earlier version of the model and better resembles the 
observed profile. 
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Figure 9. Heat Flux: Values computed from observed cloud cover, air-sea tem¬ 
perature and humidity, and solar elevation angle. 
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figure 10. Difference in SST: Values represent estimates by enhanced version 
minus estimates by earlier version. 
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Figure 11. Thermocline representation: Solid line: enhanced version estimate; 
dotted line: earlier version estimate; asterisks: observed BT data. 
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V. SUMMARY AND REMARKS 


This investigation has introduced and tested the most recent addition to the NPS 
mixed layer model, the inclusion of entrainment shear production of TKE in a finite 
entrainment zone. The study has used the NPS model to investigate the additional 
process in three ways. The first investigation was whether or not the process would 
improve annual simulations of SST. The answer to this question is not clear. In addi¬ 
tion to difficulties such as lack of data precision and the inability to estimate the possible 
effects of advection, too little is known about appropriate values for the model constants 
used to parameterize physical processes. Tuning the model to optimize performance on 
an annual time scale may inadvertently over-estimate the importance of any one process 
in order to compensate for errors in another. Since it was demonstrated that the model 
constants used in this study did not provide universally good results in tests with data 
from periods other than 1966, an optimal set of model constants has yet to be found. 
Perhaps using a model that has been tuned to better represent the structure and evolu¬ 
tion of the thermocline. rather than simply reduce SST error, would show that 
entrainment shear production provides a significant contribution to annual simulations. 
Additionally, it was noted that the apparent insignificance of entrainment shear pro¬ 
duction for annual simulations in the mid-latitudes is not indicative of its probable im¬ 
portance in operational forecasting for much shorter periods. Evidence presented here 
and in other studies suggests that this process can provide a significant contribution to 
the dynamics of the upper ocean on shorter time scales and, or at lower latitudes where 
inertial periods are longer. 

The second component of the analysis in this study was to determine the condi¬ 
tions for which entrainment shear production should provide the greatest contribution 
to the TKE of the upper ocean. The results from this analysis show that the process can 
be expected to be significant in highly stratified oceanic regions subject to variable wind 
forcing, and the process is found to be particularly significant on time-scales shorter than 
one inertial period. However, no conclusions could be reached regarding the role of 
entrainment shear production for deep mixed layers. Additionally, this particular com¬ 
ponent of the analysis primarily used two model versions which were assigned different 
parameterizations for wind shear production. No means was available to isolate the 
contribution of wind shear production for time periods of longer than about 12 hours, 
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thus it remains to be determined if entrainment shear production plays a significant role 
in an environment of sustained wind-mixing. Only the theoretical results from this study 
have suggested that the process becomes insignificant after the first inertial period fol¬ 
lowing the initiation of a given wind-forcing event. However, this analysis has extended 
the theoretical significance of the process and provided for the first time evidence of the 
theoretical conditions set forth by de Szoekc and Rhines (1976). 

Additionally, since this investigation focused on the effects of entrainment shear 
production toward improving model skill in estimating SST values, further research is 
suggested for determining its effect on mixed-layer simulations. An investigation of that 
aspect of the contribution of entrainment shear production could be launched with his¬ 
torical data collected at Station Papa, but it is not recommended due to the poor vertical 
and marginally acceptable temporal resolution of the data. 

The third aspect of this investigation was improvement in the model representation 
of the upper ocean thermal structure provided by the additional process of entrainment 
shear production. The additional process, in conjunction with the finite entrainment 
zone, was shown to provide a representation that more closely resembles observed 
thermal features. This improvement may be practically important for Navy operational 
purposes, since it provides temperature profiles with more realistic gradients at the base 
of the mixed layer, rather than discontinuities. Temperature gradients in the ocean play 
a significant role in acoustic propagation in situations tactically significant to the Navy'. 
Several additional recommendations arc made: 

• Acquiring better estimates or parameterizations of salinity, heat, and momentum 
fluxes at the surface boundary, 

• Conducting studies which account for the effects of advection. 

• Studying the effects of planetary rotation on entrainment shear production by 
conducting research at several latitudes, and particularly in tropical regions. 

• Performing similar investigations with data collected at time intervals of no longer 
than five per cent of an inertial period and with the best feasible vertical resolution 
using modern profiling instruments. 

• Concentrating research efforts to identify the effects of entrainment shear pro¬ 
duction on time scales applicable to naval operations, weather predicts,, and 
commercial fishing interests. 

• Incorporating into three-dimensional models the treatment of entrainment shear 
production in a finite entrainment zone to improve the representation of upper 
oceanic thermal characteristics over broad regions, 
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The results of this research have shown that entrainment shear production can 
play a significant role in the dynamics of the upper ocean. Additionally, the XPS model 
appears to benefit from the additional process of entrainment shear production treated 
in a finite entrainment zone. After more extensive testing and tuning, the model will be 
shown to have improved significantly the understanding of mixed layer dynamics. 
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